Dimkid Paper #1 Analysis (7-9yo and Adults)

Setup

# load libraries
library(tidyverse)
library(psych)
library(langcog) # source: https://github.com/langcog/langcog
library(RColorBrewer)
library(plotly)
library(lubridate)
library(rms)
library(lme4)
library(brms)
library(parallel)
# clear workspace
# rm(list = ls(all = T))
graphics.off()
# set number of cores for brms
n_cores <- detectCores()

# make rounding function
round2 <- function(x) {format(round(x, 2), nsmall = 2)}
# make cleanup function
cleanup <- function(datasource, age_group) {
  if(grepl("adult", age_group)) {
    
    # set target dataset
    if(datasource == "study 1"){d <- d_raw_study1}
    if(datasource == "study 1b"){d <- d_raw_study1b}
    if(datasource == "study 1c"){d <- d_raw_study1c}
    
    # enact exclusionary criteria
    d_clean_1 <- d
    
    # recode background and demographic variables
    d_clean <- d_clean_1 %>%
      mutate( # deal with study number
        study = factor(study)) %>%
      mutate( # deal with race
        race_cat2 = factor(sub(" +$", "", ethnicity)),
        race_cat3 = factor(ifelse(grepl(" ", race_cat2) == T, "multiracial",
                                  as.character(race_cat2)))) %>%
      dplyr::select(study, subid:country_selfrep, age_group, race_cat3) %>%
      rename(race_cat = race_cat3) %>%
      mutate( # deal with religion (note: only dealing with childhood religion for now)
        religion_cat2 = factor(sub(" +$", "", religionChild)),
        religion_cat3 = factor(ifelse(grepl(" ", religion_cat2) == T, 
                                      "multireligious",
                                      as.character(religion_cat2)))) %>%
      dplyr::select(study:race_cat, religion_cat3) %>%
      rename(religion_cat = religion_cat3)
    
    # remove extraneous dfs and variables
    rm(d, d_clean_1)
  }
  
  if(grepl("child", age_group)) {
    
    # set target dataset
    if(datasource == "study 2"){d <- d_raw_study2}
    # if(datasource == "study 3"){d <- d_raw_study3}
    # if(datasource == "study 4"){d <- d_raw_study4}
    
    # recode background and demographic variables
    d_clean_2 <- d %>%
      mutate( # deal with study number
        study = factor(study),
        responseNum = ifelse(!is.na(responseNum), responseNum,
                             ifelse(response == "no", 0, 
                                    ifelse(response == "kinda", 0.5, 
                                           ifelse(response == "yes", 1, NA)))))
    # NOTE: need to reconcile race/ethnicity at some point...
    # NOTE: need to deal with gender at some point...
  
    d_clean <- d_clean_2
    
    # remove extraneous dfs and variables
    rm(d, d_clean_2)
  }
  
  # remove outliers if desired
  if(chosenOutlierHandling == "remove") {
    
    d_clean <- d_clean %>%
      gather(capacity, score, happy:pride) %>%
      group_by(character, capacity) %>%
      filter(!score %in% boxplot.stats(score, 2.5)$out) %>%
      spread(capacity, score) %>%
      arrange(character, subid)
    
  }
  
  # filter characters if desired
  if(is.element("none", chosenExclude)) {} else {
    
    d_clean <- d_clean %>%
      filter(!character %in% chosenExclude)
    
    }
    
  # filter items if desired
  if(is.element("none", chosenExcludeItem)) {} else {
    d_clean <- d_clean %>%
      dplyr::filter(!capacity %in% chosenExcludeItem)
  }
  
  # drop trials <250 ms
  d_clean <- d_clean %>%
    filter(rt >= 250 | is.na(rt))
  
  # center response variable
  if(datasource == "study 1b") {
    d_clean <- d_clean %>%
      mutate(responseNumC = responseNum - 4)
  } else {
    d_clean <- d_clean %>%
      mutate(responseNumC = responseNum - 0.5)
  }
    # rename character name variables
  if("charName" %in% names(d_clean)) {
    d_clean <- d_clean %>% rename(character = charName)
  }
  
  # cleanup
  d_clean <- d_clean %>%
    filter(!is.na(subid), !is.na(character), !is.na(capacity))
  
  # return cleaned dataset
  return(d_clean)
}
# make function for stripping dataframes for dimension reducation
makeDRDF <- function(datasource, chosenCondition) {
  
  # set target dataset
  if(datasource == "study 1"){d <- d1}
  if(datasource == "study 1b"){d <- d1b}
  if(datasource == "study 1c"){d <- d1c}
  if(datasource == "study 2"){d <- d2}
  # filter by character if specified
  if(chosenCondition %in% c("beetle", "robot")) {
    d <- d %>% filter(character == chosenCondition)
  }
  # make stripped dataframe for dimension reducation analyses
  d_strip <- d %>%
    filter(!is.na(character), !is.na(subid), !is.na(capacity), capacity != "") %>%
    mutate(subid = paste(character, subid, sep = "_")) %>%
    select(subid, capacity, responseNum) %>%
    spread(capacity, responseNum) %>%
    remove_rownames() %>%
    column_to_rownames(var = "subid")
  # return stripped dataframe
  return(d_strip)
}
# make demographics functions
demoSampleSize <- function(datasource) {
  # set target dataset
  if(datasource == "study 1"){d <- d1}
  if(datasource == "study 1b"){d <- d1b}
  if(datasource == "study 1c"){d <- d1c}
  if(datasource == "study 2"){d <- d2}
  # get distinct subids
  sample_size <- d %>% distinct(subid, character) %>% count(character) %>% data.frame()
  # add total sample size  
  sample_size <- rbind(sample_size %>% mutate(character = as.character(character)),
                       c(character = "all", n = d %>% distinct(subid) %>% count() %>% as.numeric()))
  
  # return dataframe
  return(sample_size)
}
demoDuration <- function(datasource) {
  # set target dataset
  if(datasource == "study 1"){d <- d1}
  if(datasource == "study 1b"){d <- d1b}
  if(datasource == "study 1c"){d <- d1c}
  if(datasource == "study 2"){d <- d2 %>% rename(duration = sessionDuration)}
  # get sample size per character
  duration <- d %>%
    distinct(subid, character, duration) %>%
    mutate(duration = as.numeric(duration)) %>%
    group_by(character) %>%
    summarise(min_duration = min(duration, na.rm = T),
              max_duration = max(duration, na.rm = T),
              median_duration = median(duration, na.rm = T),
              mean_duration = mean(duration, na.rm = T),
              sd_duration = sd(duration, na.rm = T))
  # add total duration
  all <- d %>%
    distinct(subid, character, duration) %>%
    mutate(duration = as.numeric(duration)) %>%
    summarise(min_duration = min(duration, na.rm = T),
              max_duration = max(duration, na.rm = T),
              median_duration = median(duration, na.rm = T),
              mean_duration = mean(duration, na.rm = T),
              sd_duration = sd(duration, na.rm = T)) %>%
    mutate(character = "all")
  
  duration <- rbind(duration, all) # not sure why full_join doesn't work    
  # return dataframe
  return(duration)
}
demoAge <- function(datasource) {
  # set target dataset
  if(datasource == "study 1"){d <- d1}
  if(datasource == "study 1b"){d <- d1b}
  if(datasource == "study 1c"){d <- d1c}
  if(datasource == "study 2"){d <- d2}
  # get sample size per character
  age <- d %>%
    distinct(subid, character, age) %>%
    mutate(age = as.numeric(age)) %>%
    group_by(character) %>%
    summarise(min_age = min(age, na.rm = T),
              max_age = max(age, na.rm = T),
              median_age = median(age, na.rm = T),
              mean_age = mean(age, na.rm = T),
              sd_age = sd(age, na.rm = T))
  # add total age
  all <- d %>%
    distinct(subid, character, age) %>%
    mutate(age = as.numeric(age)) %>%
    summarise(min_age = min(age, na.rm = T),
              max_age = max(age, na.rm = T),
              median_age = median(age, na.rm = T),
              mean_age = mean(age, na.rm = T),
              sd_age = sd(age, na.rm = T)) %>%
    mutate(character = "all")
  age <- full_join(age, all)
  # return dataframe
  return(age)
}
demoGender <- function(datasource) {
  # set target dataset
  if(datasource == "study 1"){d <- d1}
  if(datasource == "study 1b"){d <- d1b}
  if(datasource == "study 1c"){d <- d1c}
  if(datasource == "study 2"){d <- d2}
  # get gender per character and overall
  gender <- data.frame(addmargins(with(d %>% distinct(subid, character, gender), 
                                       table(character, gender)))) %>%
    filter(gender != "Sum") %>%
    rename(n = Freq)
  
  gender <- gender %>%
    mutate(character = factor(ifelse(character == "Sum",
                                     "all", as.character(character)),
                              levels = c("beetle", "robot", "all"))) %>%
    arrange(character, gender) %>%
    spread(gender, n)
  
  # return dataframe
  return(gender)
}
demoRace <- function(datasource) {
  # set target dataset
  if(datasource == "study 1"){d <- d1}
  if(datasource == "study 1b"){d <- d1b}
  if(datasource == "study 1c"){d <- d1c}
  if(datasource == "study 2"){d <- d2 %>% 
    mutate(ethnicity = as.character(ethnicity),
           white = grepl("white", ethnicity) | grepl("cauc", ethnicity) | 
             grepl("euro", ethnicity) | grepl("australia", ethnicity),
           black = grepl("black", ethnicity) | grepl("africa", ethnicity),
           latinx = grepl("hisp", ethnicity) | grepl("latin", ethnicity),
           E_asian = (grepl("east asian", ethnicity) & grepl("southeast", ethnicity) == F) | 
             grepl("chin", ethnicity) | grepl("korea", ethnicity) 
           | grepl("\\(asian\\)", ethnicity) | 
             grepl("\\/asian", ethnicity),
           SSE_asian = grepl("south or southeast", ethnicity) | 
             (grepl("india", ethnicity) & grepl("native", ethnicity) == F) |
             grepl("vietnam", ethnicity) | grepl("pakistan", ethnicity),
           nat_am = grepl("native american", ethnicity),
           middle_east = grepl("middle", ethnicity)) %>%
    mutate(race_cat = ifelse(white + black + latinx + E_asian + SSE_asian + 
                               nat_am + middle_east > 1,
                             "multiracial",
                             ifelse(white, "white",
                                    ifelse(black, "black",
                                           ifelse(latinx, "latinx",
                                                  ifelse(E_asian, "east_asian",
                                                         ifelse(SSE_asian, "south_southeast_asian",
                                                                ifelse(nat_am, "native_american",
                                                                       ifelse(middle_east, "middle_eastern",
                                                                              NA)))))))))
  }
  # get race per character and overall
  race <- data.frame(addmargins(with(d %>% distinct(subid, character, race_cat), 
                                     table(character, race_cat)))) %>%
    filter(race_cat != "Sum") %>%
    rename(n = Freq)
    race <- race %>%
      mutate(character = factor(ifelse(character == "Sum",
                                       "all", as.character(character)))) %>%
      arrange(character, race_cat) %>%
      spread(race_cat, n)
  
  # return dataframe
  return(race)
}
# plotting functions
makeFacetLabs <- function(df_plotting) {
  facet_labels <- array()
  df_plotting <- df_plotting %>% mutate(character = factor(character))
  for(i in 1:length(levels(df_plotting$character))) {
    df <- df_plotting %>% filter(character == levels(df_plotting$character)[i]) %>%
      select(character, n) %>% unique()
    facet_labels[i] <- paste0(df$character, " (n = ", df$n, ")")
  }
  names(facet_labels) <- levels(df_plotting$character)
  return(facet_labels)
}

# remove outliers?
chosenOutlierHandling <- "keep"
# chosenOutlierHandling <- "remove"
# exclude any conditions (characters)?
chosenExclude <- "none"
# chosenExclude <- c("stapler", "car", "computer")
# exclude any items (mental capacities)?
# chosenExcludeItem <- "none"
# chosenExcludeItem <- "computations"
chosenExcludeItem <- c("metal", "on_off")
# NOTE: always choose minimal residual (fm = "minres") instead of ML because of non-normality
# for EFAs, what kind of correlation?
chosenCorType <- "cor" # pearson correlation
# chosenCorType <- "poly" # polychoric correlation
# for EFAs, what kind of rotation?
chosenRotType <- "oblimin" # oblimin rotation
# chosenRotType <- "oblimin" # oblimin rotation
# chosenRotType <- "none" # no rotation
data.frame("conditionsExcluded" = chosenExclude,
           "outlierHandling" = chosenOutlierHandling,
           "EFA_correlation" = chosenCorType,
           "EFA_rotation" = chosenRotType)

back to TOC

Data preparation

# study 1 (2016-07-06, adults, 2 conditions, 3-point scale, "decide what to do" and "make plans")
d_raw_study1 <- read.csv("/Users/kweisman/Documents/Research (Stanford)/Projects/Dimkid/dimkid/data/adults/us_run-01_2016-06-05_anonymized.csv") %>%
  mutate(study = "study 1", age_group = "adults") %>% select(-X)
# study 1b (2017-07-19, adults, 2 conditions, 7-point scale, "decide what to do" and "make plans")
d_raw_study1b <- read.csv("/Users/kweisman/Documents/Research (Stanford)/Projects/Dimkid/dimkid/data/adults/us_run-02_2016-07-19_anonymized.csv") %>%
  mutate(study = "study 1b", age_group = "adults") %>% select(-X)
# study 1c (2016-12-08, adults, 2 conditions, 3-point scale, "have free will" and "have intentions")
d_raw_study1c <- read.csv("/Users/kweisman/Documents/Research (Stanford)/Projects/Dimkid/dimkid/data/adults/us_run-03_2016-12-08_anonymized.csv") %>%
  mutate(study = "study 1c", age_group = "adults") %>% select(-X)
# study 2 (June - December 2016, 7-9yo, 2 conditions, 3-point-scale, "decide what to do" and "make plans")
d_raw_study2 <- read.csv("/Users/kweisman/Documents/Research (Stanford)/Projects/Dimkid/dimkid/data/children/run-01_2017-07-24_anonymized.csv") %>%
  mutate(study = "study 2", age_group = "children_79") %>% select(-X)

# clean up datasets
d1 <- cleanup("study 1", "adults")
d1b <- cleanup("study 1", "adults")
d1c <- cleanup("study 1", "adults")
d2 <- cleanup("study 2", "children")
# tweak by hand
d2 <- d2 %>%
  filter(!is.na(age)) %>%
  filter(age >= 7, age < 10) %>%
  filter(character != "elephant")

# fast rt
d_raw_study1 %>% 
  full_join(d_raw_study2) %>%
  filter(!is.na(rt)) %>%
  mutate(rt_cat = ifelse(rt < 250, "<250ms", ">=250ms")) %>% 
  count(age_group, rt_cat) %>%
  group_by(age_group) %>%
  mutate(prop = round(n/sum(n), 2))

Joining, by = c("run", "subid", "age", "gender", "ethnicity", "trialNum", "bgColor", "capacity", "capWording", "hoverTime", "rt", "response", "responseNum", "study", "age_group")
Column `run` joining factors with different levels, coercing to character vectorColumn `subid` joining factors with different levels, coercing to character vectorColumn `gender` joining factors with different levels, coercing to character vectorColumn `ethnicity` joining factors with different levels, coercing to character vectorColumn `capWording` joining factors with different levels, coercing to character vectorColumn `response` joining factors with different levels, coercing to character vector

# skipped trials
d_raw_study1 %>% 
  full_join(d_raw_study2) %>%
  filter((age >= 7 & age < 10) | (age >= 18)) %>%
  filter(!is.na(response), !response %in% c("bail", "skip")) %>%
  count(age_group, subid) %>%
  filter(n != 40) %>%
  mutate(skipped = 40 - n) %>%
  group_by(age_group) %>%
  summarise(total_skipped = sum(skipped))

Joining, by = c("run", "subid", "age", "gender", "ethnicity", "trialNum", "bgColor", "capacity", "capWording", "hoverTime", "rt", "response", "responseNum", "study", "age_group")
Column `run` joining factors with different levels, coercing to character vectorColumn `subid` joining factors with different levels, coercing to character vectorColumn `gender` joining factors with different levels, coercing to character vectorColumn `ethnicity` joining factors with different levels, coercing to character vectorColumn `capWording` joining factors with different levels, coercing to character vectorColumn `response` joining factors with different levels, coercing to character vector

# total mising
d1 %>% 
  full_join(d2) %>%
  filter((age >= 7 & age < 10) | (age >= 18)) %>%
  filter(!is.na(response), !response %in% c("bail", "skip")) %>%
  count(age_group, subid) %>%
  filter(n != 40) %>%
  mutate(missing = 40 - n) %>%
  group_by(age_group) %>%
  summarise(total_missing = sum(missing),
            prop_missing = round(total_missing/(40*200), 3))

Joining, by = c("study", "subid", "character", "age", "gender", "ethnicity", "trialNum", "bgColor", "capacity", "capWording", "hoverTime", "rt", "response", "responseNum", "age_group", "responseNumC")
Column `study` joining factors with different levels, coercing to character vectorColumn `subid` joining factors with different levels, coercing to character vectorColumn `character` joining factors with different levels, coercing to character vectorColumn `gender` joining factors with different levels, coercing to character vectorColumn `ethnicity` joining factors with different levels, coercing to character vectorColumn `capWording` joining factors with different levels, coercing to character vectorColumn `response` joining factors with different levels, coercing to character vector

# make dataframes for s1
# d1_beetle <- makeDRDF("study 1", "beetle")
# d1_robot <- makeDRDF("study 1", "robot")
d1_all <- makeDRDF("study 1", "all")
# make dataframes for follow-up studies to s1
d1b_all <- makeDRDF("study 1b", "all")
d1c_all <- makeDRDF("study 1c", "all")
# make dataframes for study 2
# d2_beetle <- makeDRDF("study 2", "beetle")
# d2_robot <- makeDRDF("study 2", "robot")
d2_all <- makeDRDF("study 2", "all")

back to TOC

Analysis plan

For all studies we conduct exploratory factor analyses using Pearson correlations to find minimum residual solutions.

For each study, we first examine maximal unrotated and rotated solutions. To determine the maximum number of factors to extract, we use the following rule of thumb: With \(p\) observations per participant, we can extract a maximum of \(k\) factors, where \((p-k)*2 > p+k\), i.e., \(k < p/3\). Thus, with 40 mental capacity items, we can extract a maximum of 13 factors.

To determine how many factors to retain, we use the following preset retention criteria, considering the unrotated maximal solution (unless otherwise noted):

• Each factor must have an eigenvalue >1.0.
• Each factor must individually account for >5% of the total variance in the maximal model.
• After rotation, each factor must be the dominant factor (i.e., the factor with the highest factor loading) for ≥1 mental capacity item.

We then examine and interpret oblimin-rotated solutions, extracting only the number of factors that meet these criteria.

Adults

Study information:

• Design: 2 conditions (beetle, robot), 3-point response scale, “decide what to do” and “have intentions”
• Population: adults
• Date conducted: 2016-07-06 (MTurk)

Demographics

# make demographics tables
demoSampleSize("study 1")

demoDuration("study 1")

demoAge("study 1")

Joining, by = c("character", "min_age", "max_age", "median_age", "mean_age", "sd_age")
Column `character` joining factor and character vector, coercing into character vector

demoGender("study 1")

demoRace("study 1") %>%
  filter(character == "all") %>%
  gather(ethnicity, count, -character) %>%
  mutate(prop = round(count/sum(count), 2)) %>% 
  arrange(desc(prop))

# gender
with(d1 %>%
       distinct(subid, character, gender) %>%
       filter(!is.na(gender)) %>%
       mutate(character = factor(character)),
     table(character, gender)) %>%
  summary()

Number of cases in table: 200 
Number of factors: 2 
Test for independence of all factors:
    Chisq = 1.9863, df = 2, p-value = 0.3704
    Chi-squared approximation may be incorrect

# age
with(d1 %>%
       distinct(subid, character, age),
     t.test(age ~ character))

    Welch Two Sample t-test

data:  age by character
t = 0.26048, df = 192.06, p-value = 0.7948
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.625759  3.424833
sample estimates:
mean in group beetle  mean in group robot 
            33.15464             32.75510 

back to TOC

Exploratory factor analysis

Step 1: Run maximal EFA (without and with rotation)

# examine scree plot
# fa.parallel(d1_all)
# run EFA without rotation with N factors
efa_d1_all_unrotated <- fa(d1_all, 13, rotate = "none",
                           cor = chosenCorType, fm = "minres")
fa.sort(efa_d1_all_unrotated)

Factor Analysis using method =  minres
Call: fa(r = d1_all, nfactors = 13, rotate = "none", fm = "minres", 
    cor = chosenCorType)
Standardized loadings (pattern matrix) based upon correlation matrix
                 MR1   MR2   MR3   MR4   MR5   MR6   MR7   MR8   MR9  MR10  MR11
happy           0.76  0.00 -0.33 -0.08 -0.22  0.08 -0.15 -0.03 -0.08 -0.04 -0.08
joy             0.76  0.01 -0.39  0.10 -0.17 -0.01 -0.08 -0.07  0.05  0.00 -0.03
love            0.75  0.11 -0.28  0.09 -0.07 -0.09 -0.07 -0.14 -0.03 -0.08  0.00
depressed       0.74  0.04 -0.37  0.04 -0.14  0.04 -0.17 -0.04 -0.21 -0.05 -0.02
fear            0.72 -0.39  0.14  0.03 -0.18  0.07  0.09  0.12 -0.03  0.13 -0.04
safe            0.71 -0.29  0.21 -0.12 -0.03 -0.06  0.04  0.06 -0.03  0.09 -0.13
tired           0.69 -0.34  0.23  0.07 -0.06  0.05  0.06  0.08  0.07 -0.07  0.09
pleasure        0.69 -0.23 -0.07  0.15 -0.15 -0.08 -0.15  0.06  0.16 -0.06  0.21
calm            0.68 -0.17  0.01 -0.08 -0.08  0.04  0.08  0.10 -0.04  0.02 -0.01
pride           0.68  0.18 -0.42  0.08  0.04  0.04  0.01 -0.08 -0.14  0.05 -0.06
desires         0.66 -0.17  0.10 -0.02 -0.03 -0.02  0.35 -0.48  0.13 -0.14 -0.08
angry           0.65 -0.04 -0.11 -0.03  0.05  0.22 -0.02  0.04 -0.10  0.17  0.02
nauseated       0.65 -0.32  0.14  0.08 -0.16  0.05  0.01 -0.06  0.14  0.03  0.03
pain            0.63 -0.52  0.19 -0.01 -0.04  0.11  0.11  0.19  0.05 -0.01  0.09
disrespected    0.63  0.06 -0.35  0.16  0.07  0.07 -0.13 -0.03 -0.07 -0.03 -0.07
guilt           0.62  0.14 -0.41  0.21  0.43  0.14  0.04  0.02  0.08  0.00 -0.03
thoughts        0.55  0.18  0.10 -0.37 -0.01  0.04 -0.04  0.10 -0.09 -0.15  0.19
embarrassed     0.52  0.14 -0.40  0.19  0.48  0.18  0.11  0.11  0.12 -0.03  0.09
odors           0.49 -0.35  0.37  0.05  0.15 -0.09 -0.03  0.10  0.07  0.01 -0.03
beliefs         0.48  0.40 -0.16 -0.14  0.04 -0.38  0.18  0.09  0.07  0.04 -0.06
self_aware      0.46  0.18  0.22 -0.30  0.09  0.00  0.12  0.13 -0.20 -0.28 -0.21
personality     0.44  0.36 -0.19 -0.13 -0.03 -0.27  0.01  0.00  0.24  0.20 -0.06
self_restraint  0.43  0.35 -0.05 -0.15  0.04 -0.19 -0.08  0.07  0.00  0.02  0.17
goal            0.41  0.21  0.19 -0.11  0.07 -0.11  0.18 -0.08  0.02 -0.03  0.07
choices         0.37  0.34  0.36 -0.20  0.06  0.09 -0.25 -0.20  0.12 -0.05  0.17
computations   -0.33  0.82 -0.07  0.14 -0.03  0.00 -0.02  0.04  0.03  0.00  0.03
recognizing     0.10  0.76  0.12  0.13 -0.21  0.15  0.11  0.11  0.15 -0.02  0.00
hungry          0.55 -0.71  0.22 -0.03  0.08 -0.02  0.06  0.02  0.03  0.03  0.05
remembering     0.14  0.66  0.16  0.10 -0.15  0.13  0.01 -0.05  0.03 -0.15  0.04
communicating   0.11  0.62  0.18  0.14 -0.18  0.30  0.11  0.18  0.15 -0.05 -0.08
intentions      0.19  0.62  0.02 -0.16  0.01  0.07  0.26  0.00 -0.20  0.00  0.09
morality        0.31  0.50 -0.13  0.02 -0.07 -0.19  0.08  0.05 -0.04  0.00  0.18
reasoning       0.34  0.44  0.31 -0.16  0.01  0.21 -0.11  0.01  0.08  0.20 -0.21
emo_recog       0.37  0.39 -0.10 -0.10  0.01 -0.27 -0.09  0.10  0.04  0.08 -0.14
seeing          0.33  0.15  0.50  0.28  0.08 -0.03 -0.07 -0.01 -0.11 -0.01  0.07
depth           0.26  0.27  0.48  0.28  0.12 -0.09 -0.16 -0.09  0.04 -0.08 -0.10
temperature     0.30  0.19  0.46  0.40  0.05 -0.22  0.04  0.00 -0.26  0.12  0.09
conscious       0.44  0.10  0.44 -0.11  0.17 -0.11 -0.24  0.06  0.08 -0.19 -0.15
sounds          0.27  0.20  0.42  0.38 -0.06 -0.02  0.06 -0.07 -0.11  0.12 -0.07
free_will       0.31  0.30  0.32 -0.40  0.15  0.22 -0.07 -0.19 -0.05  0.27  0.10
                MR12  MR13   h2   u2 com
happy          -0.04 -0.08 0.80 0.20 1.8
joy             0.09 -0.01 0.79 0.21 1.8
love            0.08  0.07 0.71 0.29 1.6
depressed       0.06 -0.03 0.79 0.21 2.0
fear           -0.16  0.15 0.82 0.18 2.3
safe           -0.04  0.07 0.68 0.32 1.8
tired          -0.11 -0.05 0.70 0.30 2.0
pleasure        0.09  0.18 0.72 0.28 2.3
calm            0.05 -0.07 0.54 0.46 1.3
pride          -0.16  0.01 0.74 0.26 2.2
desires         0.02 -0.07 0.87 0.13 2.9
angry          -0.20 -0.09 0.58 0.42 1.8
nauseated       0.04  0.05 0.61 0.39 2.0
pain           -0.01 -0.14 0.80 0.20 2.7
disrespected   -0.02  0.02 0.59 0.41 2.0
guilt          -0.03 -0.05 0.84 0.16 3.3
thoughts        0.09 -0.05 0.57 0.43 2.9
embarrassed     0.12  0.07 0.82 0.18 4.3
odors           0.11  0.05 0.56 0.44 3.4
beliefs        -0.09 -0.04 0.64 0.36 4.0
self_aware      0.15 -0.07 0.60 0.40 5.6
personality     0.09 -0.17 0.59 0.41 5.2
self_restraint -0.04 -0.14 0.43 0.57 3.5
goal           -0.14  0.15 0.36 0.64 4.0
choices        -0.07 -0.05 0.59 0.41 6.1
computations   -0.06 -0.03 0.81 0.19 1.4
recognizing    -0.02 -0.12 0.75 0.25 1.7
hungry          0.07  0.01 0.88 0.12 2.2
remembering     0.09  0.10 0.58 0.42 1.7
communicating  -0.01 -0.03 0.65 0.35 2.6
intentions     -0.01  0.18 0.60 0.40 2.3
morality       -0.10  0.10 0.46 0.54 2.9
reasoning       0.09  0.12 0.60 0.40 5.4
emo_recog       0.14  0.11 0.46 0.54 4.3
seeing          0.09 -0.08 0.50 0.50 3.1
depth          -0.13 -0.02 0.53 0.47 4.1
temperature     0.06 -0.09 0.65 0.35 4.8
conscious      -0.21  0.05 0.63 0.37 4.6
sounds          0.08  0.01 0.48 0.52 4.0
free_will       0.08  0.00 0.65 0.35 6.8

                        MR1  MR2  MR3  MR4  MR5  MR6  MR7  MR8  MR9 MR10 MR11 MR12
SS loadings           11.03 5.51 3.15 1.29 0.87 0.84 0.64 0.56 0.51 0.46 0.43 0.38
Proportion Var         0.28 0.14 0.08 0.03 0.02 0.02 0.02 0.01 0.01 0.01 0.01 0.01
Cumulative Var         0.28 0.41 0.49 0.52 0.55 0.57 0.58 0.60 0.61 0.62 0.63 0.64
Proportion Explained   0.42 0.21 0.12 0.05 0.03 0.03 0.02 0.02 0.02 0.02 0.02 0.01
Cumulative Proportion  0.42 0.64 0.76 0.81 0.84 0.87 0.90 0.92 0.94 0.96 0.97 0.99
                      MR13
SS loadings           0.32
Proportion Var        0.01
Cumulative Var        0.65
Proportion Explained  0.01
Cumulative Proportion 1.00

Mean item complexity =  3.1
Test of the hypothesis that 13 factors are sufficient.

The degrees of freedom for the null model are  780  and the objective function was  27.45 with Chi Square of  5073.12
The degrees of freedom for the model are 338  and the objective function was  2.41 

The root mean square of the residuals (RMSR) is  0.02 
The df corrected root mean square of the residuals is  0.03 

The harmonic number of observations is  196 with the empirical chi square  93.35  with prob <  1 
The total number of observations was  200  with Likelihood Chi Square =  424.01  with prob <  0.001 

Tucker Lewis Index of factoring reliability =  0.951
RMSEA index =  0.046  and the 90 % confidence intervals are  0.024 0.046
BIC =  -1366.82
Fit based upon off diagonal values = 1
Measures of factor score adequacy             
                                                MR1  MR2  MR3  MR4  MR5  MR6  MR7
Correlation of scores with factors             0.99 0.98 0.95 0.88 0.89 0.84 0.83
Multiple R square of scores with factors       0.98 0.95 0.91 0.77 0.79 0.70 0.70
Minimum correlation of possible factor scores  0.95 0.90 0.82 0.55 0.59 0.40 0.39
                                                MR8  MR9 MR10 MR11 MR12 MR13
Correlation of scores with factors             0.85 0.78 0.75 0.73 0.74 0.71
Multiple R square of scores with factors       0.73 0.61 0.57 0.54 0.54 0.50
Minimum correlation of possible factor scores  0.45 0.23 0.14 0.08 0.08 0.00

# count factors with eigenvalues > 1 and variance explained > 5%
efa_d1_all_unrotated_nfactors <- efa_d1_all_unrotated_eigenvalues %>%
  filter(SS.loadings > 1, Proportion.Explained > 0.05) %>%
  count() %>%
  as.numeric()
efa_d1_all_unrotated_nfactors

[1] 3

efa_d1_rotatedN <- fa(d1_all, efa_d1_all_unrotated_nfactors, rotate = chosenRotType,
                           cor = chosenCorType, fm = "minres")
# check that each of these factors is the dominant factor for at least one mental capacity item
efa_d1_rotatedN_loadings <- fa.sort(loadings(efa_d1_rotatedN)[]) %>%
  data.frame() %>%
  rownames_to_column("capacity") %>%
  gather(factor, loading, -capacity) %>%
  mutate(loading_abs = abs(loading)) %>%
  group_by(capacity) %>%
  top_n(1, loading_abs) %>%
  ungroup()
efa_d1_rotatedN_loadings

# drop any factors where n < 1
efa_d1_rotatedN_loadings %>% 
  count(factor) %>% 
  filter(n > 0)

# set number of factors to extract
nfactors_d1_all <- efa_d1_rotatedN_loadings %>% 
  count(factor) %>% 
  filter(n > 0) %>%
  nrow()
nfactors_d1_all

[1] 3

Step 2: Run EFA with oblimin rotation

# run EFA with rotation with N factors
efa_d1_rotatedN <- fa(d1_all, nfactors_d1_all, 
                          rotate = chosenRotType, cor = chosenCorType, fm = "minres")
fa.sort(efa_d1_rotatedN)

Factor Analysis using method =  minres
Call: fa(r = d1_all, nfactors = nfactors_d1_all, rotate = chosenRotType, 
    fm = "minres", cor = chosenCorType)
Standardized loadings (pattern matrix) based upon correlation matrix
                 MR1   MR2   MR3   h2   u2 com
pride           0.85 -0.06 -0.04 0.68 0.32 1.0
joy             0.84  0.12 -0.09 0.73 0.27 1.1
depressed       0.81  0.10 -0.06 0.68 0.32 1.0
happy           0.78  0.17 -0.04 0.68 0.32 1.1
love            0.76  0.10  0.06 0.66 0.34 1.0
guilt           0.75 -0.03 -0.05 0.53 0.47 1.0
disrespected    0.72  0.05 -0.06 0.53 0.47 1.0
embarrassed     0.66 -0.06 -0.07 0.39 0.61 1.0
beliefs         0.53 -0.17  0.24 0.40 0.60 1.6
personality     0.52 -0.18  0.19 0.34 0.66 1.5
angry           0.50  0.27  0.09 0.43 0.57 1.6
morality        0.43 -0.32  0.28 0.36 0.64 2.6
emo_recog       0.41 -0.18  0.26 0.29 0.71 2.1
self_restraint  0.40 -0.10  0.29 0.30 0.70 2.0
hungry          0.00  0.93 -0.06 0.87 0.13 1.0
computations    0.04 -0.83  0.33 0.80 0.20 1.3
pain            0.13  0.79  0.05 0.70 0.30 1.1
tired           0.17  0.70  0.20 0.64 0.36 1.3
fear            0.26  0.70  0.11 0.68 0.32 1.3
odors          -0.08  0.69  0.25 0.50 0.50 1.3
safe            0.22  0.65  0.22 0.63 0.37 1.5
nauseated       0.23  0.62  0.13 0.54 0.46 1.4
desires         0.31  0.45  0.17 0.44 0.56 2.1
calm            0.39  0.45  0.12 0.50 0.50 2.1
pleasure        0.44  0.44  0.02 0.51 0.49 2.0
depth          -0.17  0.15  0.62 0.35 0.65 1.3
reasoning       0.07 -0.02  0.60 0.39 0.61 1.0
seeing         -0.17  0.29  0.59 0.37 0.63 1.6
choices         0.03  0.08  0.59 0.36 0.64 1.0
recognizing     0.14 -0.48  0.57 0.59 0.41 2.1
remembering     0.11 -0.37  0.56 0.48 0.52 1.8
temperature    -0.13  0.21  0.55 0.30 0.70 1.4
conscious      -0.05  0.34  0.54 0.37 0.63 1.7
communicating   0.07 -0.33  0.53 0.41 0.59 1.7
sounds         -0.13  0.18  0.53 0.27 0.73 1.3
free_will       0.02  0.06  0.49 0.25 0.75 1.0
intentions      0.25 -0.38  0.43 0.41 0.59 2.6
goal            0.15  0.12  0.41 0.25 0.75 1.4
self_aware      0.15  0.18  0.40 0.27 0.73 1.7
thoughts        0.31  0.16  0.34 0.33 0.67 2.4

                       MR1  MR2  MR3
SS loadings           7.56 6.63 5.04
Proportion Var        0.19 0.17 0.13
Cumulative Var        0.19 0.35 0.48
Proportion Explained  0.39 0.34 0.26
Cumulative Proportion 0.39 0.74 1.00

 With factor correlations of 
     MR1   MR2   MR3
MR1 1.00  0.28  0.30
MR2 0.28  1.00 -0.01
MR3 0.30 -0.01  1.00

Mean item complexity =  1.5
Test of the hypothesis that 3 factors are sufficient.

The degrees of freedom for the null model are  780  and the objective function was  27.45 with Chi Square of  5073.12
The degrees of freedom for the model are 663  and the objective function was  6.67 

The root mean square of the residuals (RMSR) is  0.05 
The df corrected root mean square of the residuals is  0.05 

The harmonic number of observations is  196 with the empirical chi square  729.79  with prob <  0.036 
The total number of observations was  200  with Likelihood Chi Square =  1219.34  with prob <  2.1e-35 

Tucker Lewis Index of factoring reliability =  0.846
RMSEA index =  0.071  and the 90 % confidence intervals are  0.059 NA
BIC =  -2293.44
Fit based upon off diagonal values = 0.98
Measures of factor score adequacy             
                                                MR1  MR2  MR3
Correlation of scores with factors             0.97 0.98 0.95
Multiple R square of scores with factors       0.95 0.96 0.89
Minimum correlation of possible factor scores  0.89 0.92 0.79

# get loadings for each factor
efa_d1_rotatedN_loadings <- loadings(efa_d1_rotatedN)[] %>%
  data.frame() %>% 
  rownames_to_column(var = "capacity")

back to TOC

Factor loadings table

data.frame(loadings(fa.sort(efa_d1_rotatedN))[]) %>%
  rownames_to_column("capacity") %>%
  mutate_at(vars(starts_with("M")), funs(round2))

Children (7-9y)

Study information:

• Design: 2 conditions (beetle, robot), 3-point response scale, “decide what to do” and “have intentions”
• Population: 7-9yo children
• Date conducted: June - December 2016

Demographics

# make demographics tables
demoSampleSize("study 2")

demoDuration("study 2")

demoAge("study 2")

Joining, by = c("character", "min_age", "max_age", "median_age", "mean_age", "sd_age")
Column `character` joining factor and character vector, coercing into character vector

demoGender("study 2")

demoRace("study 2") %>%
  filter(character == "all") %>%
  gather(ethnicity, count, -character) %>%
  mutate(prop = round(count/sum(count), 2)) %>% 
  arrange(desc(prop))

# gender
with(d2 %>%
       distinct(subid, character, gender) %>%
       filter(!is.na(gender)) %>%
       mutate(character = factor(character)),
     table(character, gender)) %>%
  summary()

Number of cases in table: 197 
Number of factors: 2 
Test for independence of all factors:
    Chisq = 0.016177, df = 1, p-value = 0.8988

# age
with(d2 %>%
       distinct(subid, character, age),
     t.test(age ~ character))

    Welch Two Sample t-test

data:  age by character
t = 0.48254, df = 197.76, p-value = 0.63
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1769700  0.2916335
sample estimates:
mean in group beetle  mean in group robot 
            8.381464             8.324132 

back to TOC

Exploratory factor analysis

Step 1: Run maximal EFA (without and with rotation)

# examine scree plot
# fa.parallel(d2_all)
# run EFA without rotation with N factors
efa_d2_all_unrotated <- fa(d2_all, 13, rotate = "none",
                           cor = chosenCorType, fm = "minres")
fa.sort(efa_d2_all_unrotated)

Factor Analysis using method =  minres
Call: fa(r = d2_all, nfactors = 13, rotate = "none", fm = "minres", 
    cor = chosenCorType)
Standardized loadings (pattern matrix) based upon correlation matrix
                 MR1   MR2   MR3   MR4   MR5   MR6   MR7   MR8   MR9  MR10  MR11
happy           0.72  0.09 -0.20  0.18 -0.20 -0.06 -0.01 -0.13  0.09 -0.03  0.12
joy             0.70  0.00 -0.26  0.10 -0.18 -0.10 -0.03 -0.05  0.06 -0.01  0.05
depressed       0.69 -0.11 -0.13  0.05 -0.06 -0.13  0.12  0.06  0.21 -0.09 -0.15
disrespected    0.69  0.00 -0.11 -0.07 -0.03 -0.02  0.16 -0.16 -0.09  0.03  0.10
pride           0.68  0.14 -0.31  0.06 -0.05 -0.02  0.01 -0.04 -0.01 -0.04 -0.03
love            0.60  0.03 -0.20  0.11  0.03 -0.11  0.09 -0.19 -0.05  0.05  0.07
pleasure        0.60  0.04 -0.18  0.10 -0.23 -0.17 -0.31  0.04  0.01 -0.04  0.08
safe            0.58 -0.08  0.23  0.06 -0.06 -0.06 -0.19  0.16 -0.24  0.00 -0.10
guilt           0.58  0.07 -0.18  0.13  0.25  0.19 -0.07 -0.10 -0.07 -0.08  0.14
thoughts        0.57  0.01  0.18 -0.05 -0.11 -0.16 -0.03  0.03 -0.09  0.03 -0.03
angry           0.57 -0.17 -0.04 -0.07  0.05  0.16  0.11  0.06 -0.04  0.01 -0.04
fear            0.55 -0.37  0.11  0.11 -0.08  0.03  0.01  0.17  0.09  0.03 -0.06
calm            0.55 -0.07  0.01  0.02 -0.10 -0.09 -0.14  0.05 -0.01 -0.18  0.10
personality     0.54  0.29 -0.01 -0.12  0.20  0.12  0.01  0.18 -0.01  0.06 -0.09
beliefs         0.54  0.23 -0.07 -0.09  0.05  0.24 -0.12 -0.05 -0.01  0.13  0.06
embarrassed     0.53  0.03 -0.30  0.10  0.22  0.34  0.23 -0.06 -0.01 -0.03  0.02
desires         0.52 -0.19 -0.03 -0.01 -0.08  0.21 -0.07 -0.04  0.02  0.10  0.03
free_will       0.49  0.00  0.32 -0.22 -0.09  0.11 -0.21 -0.18  0.08  0.03 -0.02
morality        0.44  0.37 -0.13 -0.08  0.21 -0.26 -0.05  0.26 -0.18  0.33  0.07
choices         0.43  0.05  0.26 -0.38 -0.07  0.20 -0.05 -0.13  0.11  0.05 -0.15
tired           0.39 -0.35  0.07  0.11  0.16 -0.16  0.24  0.06  0.12 -0.07 -0.11
goal            0.35  0.31 -0.03 -0.09 -0.06  0.10 -0.11  0.22  0.02 -0.18 -0.25
computations   -0.01  0.80 -0.04  0.02  0.08  0.10 -0.09  0.03  0.10 -0.11  0.08
hungry          0.38 -0.77  0.22  0.02  0.14 -0.04  0.01 -0.08  0.06  0.05  0.00
pain            0.45 -0.64  0.21 -0.06  0.07 -0.09 -0.05 -0.04 -0.01  0.10 -0.01
remembering     0.06  0.58  0.17  0.13 -0.05  0.00  0.01 -0.17  0.19  0.28 -0.24
odors           0.15 -0.53  0.38  0.06  0.09  0.12 -0.06  0.16  0.06  0.03  0.19
emo_recog       0.34  0.50 -0.05  0.09  0.11 -0.07  0.22  0.14 -0.09  0.07 -0.02
nauseated       0.30 -0.43  0.06  0.11  0.27  0.05  0.00  0.07  0.16  0.06 -0.10
intentions      0.35  0.40  0.15 -0.24  0.11  0.03  0.02 -0.12  0.09 -0.14  0.08
recognizing     0.20  0.32  0.13  0.10  0.11  0.00 -0.11  0.21  0.06  0.00  0.04
communicating   0.09  0.30  0.16  0.30 -0.17  0.03  0.19  0.11  0.29  0.04  0.04
conscious       0.36  0.08  0.48  0.21 -0.14  0.01  0.16 -0.12 -0.27 -0.14 -0.04
self_aware      0.27  0.21  0.46 -0.03 -0.10  0.09  0.23  0.00 -0.31 -0.10 -0.09
temperature    -0.05  0.34  0.41  0.31  0.34 -0.26 -0.15 -0.23  0.08 -0.15 -0.09
sounds         -0.07  0.10  0.40  0.10 -0.18  0.23  0.08  0.23  0.14 -0.08  0.25
depth           0.14  0.22  0.36  0.12  0.34  0.02 -0.14 -0.07 -0.02  0.02  0.17
reasoning       0.23  0.28  0.36  0.05 -0.16 -0.01 -0.05  0.01  0.05  0.12 -0.06
self_restraint  0.34  0.19  0.15 -0.56  0.11 -0.33  0.20  0.05  0.17 -0.16  0.15
seeing         -0.07  0.13  0.26 -0.03 -0.19 -0.08  0.19 -0.08  0.02  0.27  0.20
                MR12  MR13   h2   u2 com
happy           0.05 -0.09 0.69 0.31 1.8
joy             0.11  0.02 0.62 0.38 1.6
depressed      -0.01  0.00 0.63 0.37 1.7
disrespected    0.09  0.08 0.58 0.42 1.4
pride           0.01 -0.15 0.61 0.39 1.7
love            0.03  0.14 0.50 0.50 1.9
pleasure       -0.09 -0.02 0.60 0.40 2.5
safe            0.09 -0.22 0.59 0.41 2.8
guilt          -0.05 -0.03 0.53 0.47 2.4
thoughts        0.15  0.08 0.44 0.56 1.8
angry           0.02 -0.16 0.44 0.56 1.7
fear           -0.25 -0.09 0.58 0.42 2.9
calm           -0.02  0.10 0.40 0.60 1.8
personality     0.00  0.09 0.50 0.50 2.6
beliefs         0.08  0.21 0.50 0.50 2.8
embarrassed    -0.13 -0.05 0.63 0.37 3.7
desires        -0.17  0.15 0.42 0.58 2.3
free_will      -0.03 -0.10 0.51 0.49 3.5
morality        0.03 -0.09 0.69 0.31 5.9
choices         0.16 -0.09 0.54 0.46 4.6
tired           0.12  0.05 0.45 0.55 4.7
goal           -0.04  0.14 0.42 0.58 5.4
computations    0.10  0.00 0.71 0.29 1.2
hungry          0.03  0.01 0.81 0.19 1.8
pain            0.06  0.08 0.70 0.30 2.3
remembering    -0.03 -0.08 0.59 0.41 2.8
odors           0.09 -0.06 0.55 0.45 3.0
emo_recog       0.13  0.05 0.49 0.51 3.1
nauseated       0.02  0.08 0.42 0.58 3.6
intentions      0.05 -0.12 0.43 0.57 4.2
recognizing    -0.14  0.05 0.26 0.74 4.7
communicating  -0.02 -0.12 0.40 0.60 5.9
conscious      -0.12  0.01 0.58 0.42 4.2
self_aware     -0.03  0.01 0.52 0.48 4.2
temperature     0.08  0.04 0.68 0.32 6.4
sounds          0.23  0.06 0.47 0.53 5.6
depth          -0.15 -0.07 0.41 0.59 4.8
reasoning      -0.11  0.22 0.37 0.63 4.8
self_restraint -0.19  0.01 0.76 0.24 4.5
seeing         -0.08  0.02 0.30 0.70 6.0

                       MR1  MR2  MR3  MR4  MR5  MR6  MR7  MR8  MR9 MR10 MR11 MR12
SS loadings           8.25 4.13 2.18 1.05 0.96 0.84 0.73 0.66 0.62 0.54 0.52 0.44
Proportion Var        0.21 0.10 0.05 0.03 0.02 0.02 0.02 0.02 0.02 0.01 0.01 0.01
Cumulative Var        0.21 0.31 0.36 0.39 0.41 0.44 0.45 0.47 0.49 0.50 0.51 0.52
Proportion Explained  0.39 0.19 0.10 0.05 0.05 0.04 0.03 0.03 0.03 0.03 0.02 0.02
Cumulative Proportion 0.39 0.58 0.68 0.73 0.78 0.82 0.85 0.88 0.91 0.94 0.96 0.98
                      MR13
SS loadings           0.40
Proportion Var        0.01
Cumulative Var        0.53
Proportion Explained  0.02
Cumulative Proportion 1.00

Mean item complexity =  3.4
Test of the hypothesis that 13 factors are sufficient.

The degrees of freedom for the null model are  780  and the objective function was  17.08 with Chi Square of  3157.08
The degrees of freedom for the model are 338  and the objective function was  1.79 

The root mean square of the residuals (RMSR) is  0.02 
The df corrected root mean square of the residuals is  0.03 

The harmonic number of observations is  198 with the empirical chi square  142.71  with prob <  1 
The total number of observations was  200  with Likelihood Chi Square =  314.7  with prob <  0.81 

Tucker Lewis Index of factoring reliability =  1.024
RMSEA index =  0.016  and the 90 % confidence intervals are  0 0.018
BIC =  -1476.13
Fit based upon off diagonal values = 0.99
Measures of factor score adequacy             
                                                MR1  MR2  MR3  MR4  MR5  MR6  MR7
Correlation of scores with factors             0.98 0.96 0.91 0.86 0.83 0.83 0.79
Multiple R square of scores with factors       0.95 0.93 0.83 0.74 0.69 0.68 0.62
Minimum correlation of possible factor scores  0.90 0.86 0.67 0.48 0.38 0.36 0.23
                                                MR8  MR9 MR10 MR11 MR12  MR13
Correlation of scores with factors             0.77 0.77 0.75 0.73 0.71  0.67
Multiple R square of scores with factors       0.60 0.59 0.56 0.53 0.50  0.45
Minimum correlation of possible factor scores  0.19 0.18 0.12 0.06 0.00 -0.10

# count factors with eigenvalues > 1 and variance explained > 5%
efa_d2_all_unrotated_nfactors <- efa_d2_all_unrotated_eigenvalues %>%
  filter(SS.loadings > 1, Proportion.Explained > 0.05) %>%
  count() %>%
  as.numeric()
efa_d2_all_unrotated_nfactors

[1] 3

efa_d2_rotatedN <- fa(d2_all, efa_d2_all_unrotated_nfactors, rotate = chosenRotType,
                           cor = chosenCorType, fm = "minres")
# check that each of these factors is the dominant factor for at least one mental capacity item
efa_d2_rotatedN_loadings <- fa.sort(loadings(efa_d2_rotatedN)[]) %>%
  data.frame() %>%
  rownames_to_column("capacity") %>%
  gather(factor, loading, -capacity) %>%
  mutate(loading_abs = abs(loading)) %>%
  group_by(capacity) %>%
  top_n(1, loading_abs) %>%
  ungroup()
efa_d2_rotatedN_loadings

# drop any factors where n < 1
efa_d2_rotatedN_loadings %>% 
  count(factor) %>% 
  filter(n > 0)

# set number of factors to extract
nfactors_d2_all <- efa_d2_rotatedN_loadings %>% 
  count(factor) %>% 
  filter(n > 0) %>%
  nrow()
nfactors_d2_all

[1] 3

Step 2: Run EFA with oblimin rotation

# run EFA with rotation with N factors
efa_d2_rotatedN <- fa(d2_all, nfactors_d2_all, 
                          rotate = chosenRotType, cor = chosenCorType, fm = "minres")
fa.sort(efa_d2_rotatedN)

Factor Analysis using method =  minres
Call: fa(r = d2_all, nfactors = nfactors_d2_all, rotate = chosenRotType, 
    fm = "minres", cor = chosenCorType)
Standardized loadings (pattern matrix) based upon correlation matrix
                 MR1   MR2   MR3    h2   u2 com
pride           0.81 -0.11 -0.08 0.594 0.41 1.1
joy             0.76  0.05 -0.07 0.555 0.44 1.0
happy           0.74  0.00  0.02 0.555 0.44 1.0
disrespected    0.65  0.11  0.08 0.488 0.51 1.1
depressed       0.65  0.20  0.02 0.504 0.50 1.2
love            0.64  0.03 -0.03 0.403 0.60 1.0
embarrassed     0.62 -0.02 -0.12 0.346 0.65 1.1
pleasure        0.61  0.03  0.00 0.375 0.63 1.0
guilt           0.60 -0.01  0.01 0.358 0.64 1.0
beliefs         0.51 -0.12  0.16 0.336 0.66 1.3
personality     0.49 -0.15  0.24 0.372 0.63 1.7
angry           0.48  0.26  0.05 0.356 0.64 1.6
morality        0.47 -0.26  0.13 0.306 0.69 1.8
calm            0.44  0.17  0.12 0.299 0.70 1.5
desires         0.42  0.27  0.04 0.301 0.70 1.7
thoughts        0.35  0.17  0.34 0.364 0.64 2.5
goal            0.35 -0.21  0.17 0.212 0.79 2.2
hungry          0.06  0.87  0.04 0.786 0.21 1.0
pain            0.15  0.77  0.10 0.662 0.34 1.1
computations    0.12 -0.76  0.25 0.650 0.35 1.3
odors          -0.20  0.64  0.22 0.425 0.57 1.4
fear            0.33  0.49  0.12 0.441 0.56 1.9
nauseated       0.15  0.47 -0.02 0.270 0.73 1.2
remembering     0.02 -0.44  0.36 0.333 0.67 1.9
tired           0.23  0.42  0.04 0.266 0.73 1.6
emo_recog       0.38 -0.40  0.22 0.360 0.64 2.6
self_aware     -0.05  0.02  0.57 0.303 0.70 1.0
conscious       0.00  0.16  0.56 0.331 0.67 1.2
reasoning      -0.02 -0.08  0.51 0.265 0.73 1.0
depth          -0.09 -0.04  0.45 0.182 0.82 1.1
free_will       0.19  0.21  0.43 0.331 0.67 1.9
temperature    -0.23 -0.17  0.43 0.212 0.79 1.9
sounds         -0.28  0.03  0.37 0.144 0.86 1.9
intentions      0.24 -0.24  0.37 0.293 0.71 2.5
choices         0.19  0.13  0.36 0.231 0.77 1.8
safe            0.31  0.27  0.34 0.380 0.62 2.9
recognizing     0.12 -0.20  0.29 0.154 0.85 2.1
seeing         -0.21 -0.04  0.27 0.083 0.92 1.9
communicating   0.01 -0.20  0.26 0.111 0.89 1.9
self_restraint  0.21 -0.04  0.26 0.144 0.86 2.0

                       MR1  MR2  MR3
SS loadings           6.84 4.15 3.09
Proportion Var        0.17 0.10 0.08
Cumulative Var        0.17 0.27 0.35
Proportion Explained  0.49 0.29 0.22
Cumulative Proportion 0.49 0.78 1.00

 With factor correlations of 
     MR1   MR2   MR3
MR1 1.00  0.16  0.34
MR2 0.16  1.00 -0.02
MR3 0.34 -0.02  1.00

Mean item complexity =  1.6
Test of the hypothesis that 3 factors are sufficient.

The degrees of freedom for the null model are  780  and the objective function was  17.08 with Chi Square of  3157.08
The degrees of freedom for the model are 663  and the objective function was  4.79 

The root mean square of the residuals (RMSR) is  0.05 
The df corrected root mean square of the residuals is  0.05 

The harmonic number of observations is  198 with the empirical chi square  788.89  with prob <  0.00052 
The total number of observations was  200  with Likelihood Chi Square =  874.92  with prob <  5.4e-08 

Tucker Lewis Index of factoring reliability =  0.894
RMSEA index =  0.047  and the 90 % confidence intervals are  0.032 0.047
BIC =  -2637.86
Fit based upon off diagonal values = 0.95
Measures of factor score adequacy             
                                                MR1  MR2  MR3
Correlation of scores with factors             0.96 0.96 0.90
Multiple R square of scores with factors       0.92 0.92 0.81
Minimum correlation of possible factor scores  0.85 0.84 0.62

# get loadings for each factor
efa_d2_rotatedN_loadings <- loadings(efa_d2_rotatedN)[] %>%
  data.frame() %>% 
  rownames_to_column(var = "capacity")

back to TOC

Factor loadings table

data.frame(loadings(fa.sort(efa_d2_rotatedN))[]) %>%
  rownames_to_column("capacity") %>%
  mutate_at(vars(starts_with("M")), funs(round2))

Big factor loadings table

# manually set 3 factors
order_s1_manual <- loadings(fa.sort(fa(d1_all, nfactors = 3,
                                       rotate = chosenRotType, cor = chosenCorType)))[] %>%
  data.frame() %>%
  rownames_to_column(var = "capacity") %>%
  rownames_to_column(var = "order1_manual") %>%
  rename(s1_heart = MR1, s1_body = MR2, s1_mind = MR3)
order_s1 <- loadings(fa.sort(efa_d1_rotatedN))[] %>%
  data.frame() %>%
  rownames_to_column(var = "capacity") %>%
  rownames_to_column(var = "order1") %>%
  rename(s1_heart = MR1, s1_body = MR2, s1_mind = MR3)
order_s2 <- loadings(fa.sort(efa_d2_rotatedN))[] %>%
  data.frame() %>%
  rownames_to_column(var = "capacity") %>%
  rename(s2_body = MR2, s2_heart = MR1, s2_mind = MR3)
bigTable <- order_s1 %>% 
  # order_s1_manual %>% # could substitute
  full_join(order_s2) %>%
  mutate_at(vars(starts_with("s")), funs(round2)) %>%
  select(order1, # could subistitute order_s1
  # select(order1_manual, # could subistitute order_s1
         capacity, ends_with("heart"), ends_with("body"), ends_with("mind")) 

Joining, by = "capacity"

bigTable

Figures

3D scatterplots

Factor loadings for the 40 mental capacities on the three rotated factors in Study 1. Items are colored by their dominant factor loading: Items that loaded most strongly on the body factor (bodily states and will) are in red; items that loaded most strongly on the heart factor (social-emotional experiences and morality) are in blue; and items that loaded most strongly on the mind factor (perceptual-cognitive abilities and goal pursuit) are in green.

Study 1 (adults)

# set up labels for plot (shortened version of mental capacity items)
wording_s1 <- loadings(efa_d1_rotatedN)[] %>%
  data.frame() %>%
  rownames_to_column(var = "item") %>%
  select(item) %>%
  mutate(wording = factor(
    recode(item,
           happy = "feel happy",
           depressed = "feel sad",
           fear = "feel scared",
           angry = "get angry",
           calm = "feel calm",
           sounds = "hear sounds",
           seeing = "see things",
           temperature = "sense temperatures",
           odors = "smell things",
           depth = "sense whether something is close by or far away",
           computations = "do math",
           thoughts = "have thoughts",
           reasoning = "figure out how to do things",
           remembering = "remember things",
           beliefs = "have beliefs, like when you think something is true",
           hungry = "get hungry",
           tired = "feel tired",
           pain = "feel pain",
           nauseated = "feel sick, like when you feel like you might throw up",
           safe = "feel safe",
           love = "feel love",
           recognizing = "recognize somebody else",
           communicating = "communicate with somebody else",
           guilt = "feel guilty",
           disrespected = "get hurt feelings",
           free_will = "decide what to do",
           choices = "make choices",
           self_restraint = "have self-control, like when you stop yourself from doing something you shouldn't do",
           intentions = "make plans",
           goal = "have goals, like when you're working hard to do something or make something happen",
           conscious = "be aware of things",
           self_aware = "be aware of itself",
           desires = "have desires, like when you really want something",
           embarrassed = "feel embarrassed",
           emo_recog = "understand how somebody else is feeling",
           joy = "feel joy",
           morality = "know what's nice and what's mean",
           personality = "have a personality, like when someone is shy and somebody else is silly",
           pleasure = "feel pleasure, like when something feels really good",
           pride = "feel proud"))) %>%
    mutate(short = factor(
      recode(item,
             happy = "happy",
             depressed = "sad",
           fear = "scared",
           angry = "angry",
           calm = "calm",
           sounds = "hear",
           seeing = "see",
           temperature = "temperatures",
           odors = "smell",
           depth = "depth",
           computations = "math",
           thoughts = "thoughts",
           reasoning = "figure out",
           remembering = "remember",
           beliefs = "beliefs",
           hungry = "hungry",
           tired = "tired",
           pain = "pain",
           nauseated = "sick",
           safe = "safe",
           love = "love",
           recognizing = "recognize",
           communicating = "communicate",
           guilt = "guilty",
           disrespected = "hurt feelings",
           free_will = "decide",
           choices = "choices",
           self_restraint = "self-control",
           intentions = "plans",
           goal = "goals",
           conscious = "aware",
           self_aware = "self-aware",
           desires = "desires",
           embarrassed = "embarrassed",
           emo_recog = "empathy",
           joy = "joy",
           morality = "morality",
           personality = "personality",
           pleasure = "pleasure",
           pride = "pride")))
# make dataframe for plotting
scatter_plotting <- loadings(efa_d1_rotatedN)[] %>%
  data.frame() %>%
  rownames_to_column(var = "item") %>%
  rename(HEART = MR1,
         BODY = MR2,
         MIND = MR3) %>%
  full_join(wording_s1) %>%
  mutate(dominant = factor(
    ifelse(pmax(abs(BODY), abs(HEART), abs(MIND)) == abs(BODY), "BODY",
           ifelse(pmax(abs(BODY), abs(HEART), abs(MIND)) == abs(HEART), "HEART",
                  ifelse(pmax(abs(BODY), abs(HEART), abs(MIND)) == abs(MIND), "MIND",
                         NA)))),
    size = ifelse(pmax(abs(BODY), abs(HEART), abs(MIND)) == abs(BODY), abs(BODY),
                  ifelse(pmax(abs(BODY), abs(HEART), abs(MIND)) == abs(HEART), abs(HEART),
                         ifelse(pmax(abs(BODY), abs(HEART), abs(MIND)) == abs(MIND), abs(MIND),
                                NA))),
    color = ifelse(dominant == "BODY", "#E41A1C",
                   ifelse(dominant == "HEART", "#377EB8",
                          ifelse(dominant == "MIND", "#4DAF4A",
                                 NA))))

Joining, by = "item"

# plot!
figS1 <- plot_ly(scatter_plotting, x = ~HEART, y = ~BODY, z = ~MIND,
             type = "scatter3d",
             color = ~dominant, colors = c("#E41A1C", "#377EB8", "#4DAF4A"),
             marker = list(size = 4),
             text = ~short,
             textfont = list(size = 15),
             mode = "text+markers",
             showlegend = TRUE)
figS1

Study 2 (7-9yo)

# set up labels for plot (shortened version of mental capacity items)
wording_s2 <- loadings(efa_d2_rotatedN)[] %>%
  data.frame() %>%
  rownames_to_column(var = "item") %>%
  select(item) %>%
  mutate(wording = factor(
    recode(item,
           happy = "feel happy",
           depressed = "feel sad",
           fear = "feel scared",
           angry = "get angry",
           calm = "feel calm",
           sounds = "hear sounds",
           seeing = "see things",
           temperature = "sense temperatures",
           odors = "smell things",
           depth = "sense whether something is close by or far away",
           computations = "do math",
           thoughts = "have thoughts",
           reasoning = "figure out how to do things",
           remembering = "remember things",
           beliefs = "have beliefs, like when you think something is true",
           hungry = "get hungry",
           tired = "feel tired",
           pain = "feel pain",
           nauseated = "feel sick, like when you feel like you might throw up",
           safe = "feel safe",
           love = "feel love",
           recognizing = "recognize somebody else",
           communicating = "communicate with somebody else",
           guilt = "feel guilty",
           disrespected = "get hurt feelings",
           free_will = "decide what to do",
           choices = "make choices",
           self_restraint = "have self-control, like when you stop yourself from doing something you shouldn't do",
           intentions = "make plans",
           goal = "have goals, like when you're working hard to do something or make something happen",
           conscious = "be aware of things",
           self_aware = "be aware of itself",
           desires = "have desires, like when you really want something",
           embarrassed = "feel embarrassed",
           emo_recog = "understand how somebody else is feeling",
           joy = "feel joy",
           morality = "know what's nice and what's mean",
           personality = "have a personality, like when someone is shy and somebody else is silly",
           pleasure = "feel pleasure, like when something feels really good",
           pride = "feel proud"))) %>%
    mutate(short = factor(
      recode(item,
             happy = "happy",
             depressed = "sad",
           fear = "scared",
           angry = "angry",
           calm = "calm",
           sounds = "hear",
           seeing = "see",
           temperature = "temperatures",
           odors = "smell",
           depth = "depth",
           computations = "math",
           thoughts = "thoughts",
           reasoning = "figure out",
           remembering = "remember",
           beliefs = "beliefs",
           hungry = "hungry",
           tired = "tired",
           pain = "pain",
           nauseated = "sick",
           safe = "safe",
           love = "love",
           recognizing = "recognize",
           communicating = "communicate",
           guilt = "guilty",
           disrespected = "hurt feelings",
           free_will = "decide",
           choices = "choices",
           self_restraint = "self-control",
           intentions = "plans",
           goal = "goals",
           conscious = "aware",
           self_aware = "self-aware",
           desires = "desires",
           embarrassed = "embarrassed",
           emo_recog = "empathy",
           joy = "joy",
           morality = "morality",
           personality = "personality",
           pleasure = "pleasure",
           pride = "pride")))
# make dataframe for plotting
scatter_plotting <- loadings(efa_d2_rotatedN)[] %>%
  data.frame() %>%
  rownames_to_column(var = "item") %>%
  rename(HEART = MR1,
         BODY = MR2,
         MIND = MR3) %>%
  full_join(wording_s2) %>%
  mutate(dominant = factor(
    ifelse(pmax(abs(BODY), abs(HEART), abs(MIND)) == abs(BODY), "BODY",
           ifelse(pmax(abs(BODY), abs(HEART), abs(MIND)) == abs(HEART), "HEART",
                  ifelse(pmax(abs(BODY), abs(HEART), abs(MIND)) == abs(MIND), "MIND",
                         NA)))),
    size = ifelse(pmax(abs(BODY), abs(HEART), abs(MIND)) == abs(BODY), abs(BODY),
                  ifelse(pmax(abs(BODY), abs(HEART), abs(MIND)) == abs(HEART), abs(HEART),
                         ifelse(pmax(abs(BODY), abs(HEART), abs(MIND)) == abs(MIND), abs(MIND),
                                NA))),
    color = ifelse(dominant == "BODY", "#E41A1C",
                   ifelse(dominant == "HEART", "#4DAF4A",
                          ifelse(dominant == "MIND", "#E41A1C",
                                 NA))))

Joining, by = "item"

# plot!
figS2 <- plot_ly(scatter_plotting, x = ~HEART, y = ~BODY, z = ~MIND,
             type = "scatter3d",
             color = ~dominant, colors = c("#E41A1C", "#377EB8", "#4DAF4A"),
             marker = list(size = 4),
             text = ~short,
             textfont = list(size = 15),
             mode = "text+markers",
             showlegend = TRUE)
figS2

Heatmaps

NOTE: set to 3 factors manually, for now.

Study 1 (adults)

factors_s1 <- fa.sort(fa(d1_all, nfactors = 3, cor = chosenCorType, rotate = chosenRotType)$loadings[]) %>%
  data.frame() %>%
  rownames_to_column(var = "item") %>%
  full_join(wording_s1) %>%
  select(wording, MR1, MR2, MR3) %>%
  rename(capacity = wording, Factor1 = MR1, Factor2 = MR2, Factor3 = MR3) %>%
  rownames_to_column(var = "order") %>%
  mutate(order = as.numeric(order))

Joining, by = "item"

factors_s1_long <- factors_s1 %>%
  gather(factor, loading, -capacity, -order) %>%
  mutate(factor = factor(gsub("Factor", "F", factor))) %>%
  # mutate(factor = factor(gsub("Factor", "F", factor),
  #                        levels = c("F1", "F3", "F2"))) %>%
  # mutate(factor = factor(gsub("Factor", "F", factor),
  #                        levels = c("F2", "F1", "F3"))) %>%
  arrange(order, factor)
factors_s1_blank1 <- factors_s1_long %>%
  mutate(loading = rep(100, length(factors_s1_long$loading)))
# factors_s1_blank2 <- factors_s1_long %>%
#   mutate(loading = ifelse(factor == "F1", loading, rep(100, length(factors_s1_long$loading)*2/3)))
factors_s1_blank2 <- factors_s1_long %>%
  mutate(loading = ifelse(factor == "F2", loading, rep(100, length(factors_s1_long$loading)*2/3)))
factors_s1_blank3 <- factors_s1_long %>%
  mutate(loading = ifelse(factor != "F3", loading, rep(100, length(factors_s1_long$loading)*1/3)))
# ggplot(factors_s1_blank1, aes(x = factor,
#                             y = reorder(capacity, desc(order)), fill = loading)) +
#   geom_tile(color = "black") +
#   # geom_text(aes(label = format(round(loading, 2), nsmall = 2))) +
#   scale_fill_distiller(palette = "RdYlBu", limits = c(-1, 1), breaks = c(-1, 0, 1),
#                        guide = guide_colorbar(title = element_blank(),
#                                               barheight = 20)) +
#   scale_x_discrete(position = "top") +
#   theme_minimal() +
#   theme(text = element_text(size = 24),
#         axis.text.x = element_text(size = 28),
#         axis.title = element_blank(),
#         panel.grid = element_blank()) # 1000 by 1000
# 
# ggplot(factors_s1_blank2, aes(x = factor,
#                               y = reorder(capacity, desc(order)), fill = loading)) +
#   geom_tile(color = "black") +
#   # geom_text(aes(label = format(round(loading, 2), nsmall = 2))) +
#   scale_fill_distiller(palette = "RdYlBu", limits = c(-1, 1), breaks = c(-1, 0, 1),
#                        guide = guide_colorbar(title = element_blank(),
#                                               barheight = 20)) +
#   scale_x_discrete(position = "top") +
#   theme_minimal() +
#   theme(text = element_text(size = 24),
#         axis.text.x = element_text(size = 28),
#         axis.title = element_blank(),
#         panel.grid = element_blank()) # 1000 by 1000
# 
# ggplot(factors_s1_blank3, aes(x = factor,
#                               y = reorder(capacity, desc(order)), fill = loading)) +
#   geom_tile(color = "black") +
#   # geom_text(aes(label = format(round(loading, 2), nsmall = 2))) +
#   scale_fill_distiller(palette = "RdYlBu", limits = c(-1, 1), breaks = c(-1, 0, 1),
#                        guide = guide_colorbar(title = element_blank(),
#                                               barheight = 20)) +
#   scale_x_discrete(position = "top") +
#   theme_minimal() +
#   theme(text = element_text(size = 24),
#         axis.text.x = element_text(size = 28),
#         axis.title = element_blank(),
#         panel.grid = element_blank()) # 1000 by 1000
ggplot(factors_s1_long, aes(x = factor,
                            y = reorder(capacity, desc(order)), fill = loading)) +
  geom_tile(color = "black") +
  geom_text(aes(label = format(round(loading, 2), nsmall = 2)), size = 6) +
  scale_fill_distiller(palette = "RdYlBu", limits = c(-1, 1), breaks = c(-1, 0, 1),
                       guide = guide_colorbar(title = element_blank(),
                                              barheight = 20)) +
  scale_x_discrete(position = "top") +
  # geom_rect(aes(xmin = 0.51, xmax = 1.49, ymin = 14.55, ymax = 20.45),
  #           alpha = 0, color = "black", size = .5) +
  # geom_rect(aes(xmin = 1.51, xmax = 2.49, ymin = 6.55, ymax = 14.45),
  #           alpha = 0, color = "black", size = .5) +
  # geom_rect(aes(xmin = 2.51, xmax = 3.49, ymin = 0.55, ymax = 6.45),
  #           alpha = 0, color = "black", size = .5) +
  # theme_bw() +
  theme_minimal() +
  theme(text = element_text(size = 24),
        axis.text.x = element_text(size = 28),
        axis.title = element_blank(),
        panel.grid = element_blank()) # 1000 by 1000

Study 2 (7-9y)

factors_s2 <- fa.sort(fa(d2_all, nfactors = 3, cor = chosenCorType, rotate = chosenRotType)$loadings[]) %>%
  data.frame() %>%
  rownames_to_column(var = "item") %>%
  full_join(wording_s2) %>%
  select(wording, MR1, MR2, MR3) %>%
  rename(capacity = wording, Factor1 = MR1, Factor2 = MR2, Factor3 = MR3) %>%
  rownames_to_column(var = "order") %>%
  mutate(order = as.numeric(order))

Joining, by = "item"

factors_s2_long <- factors_s2 %>%
  gather(factor, loading, -capacity, -order) %>%
  mutate(factor = factor(gsub("Factor", "F", factor))) %>%
  # mutate(factor = factor(gsub("Factor", "F", factor),
  #                        levels = c("F1", "F3", "F2"))) %>%
  # mutate(factor = factor(gsub("Factor", "F", factor),
  #                        levels = c("F2", "F1", "F3"))) %>%
  arrange(order, factor)
factors_s2_blank1 <- factors_s2_long %>%
  mutate(loading = rep(100, length(factors_s2_long$loading)))
factors_s2_blank2 <- factors_s2_long %>%
  mutate(loading = ifelse(factor == "F1", loading, rep(100, length(factors_s2_long$loading)*2/3)))
# factors_s2_blank2 <- factors_s2_long %>%
#   mutate(loading = ifelse(factor == "F2", loading, rep(100, length(factors_s2_long$loading)*2/3)))
factors_s2_blank3 <- factors_s2_long %>%
  mutate(loading = ifelse(factor != "F3", loading, rep(100, length(factors_s2_long$loading)*1/3)))
# ggplot(factors_s2_blank1, aes(x = factor,
#                             y = reorder(capacity, desc(order)), fill = loading)) +
#   geom_tile(color = "black") +
#   # geom_text(aes(label = format(round(loading, 2), nsmall = 2))) +
#   scale_fill_distiller(palette = "RdYlBu", limits = c(-1, 1), breaks = c(-1, 0, 1),
#                        guide = guide_colorbar(title = element_blank(),
#                                               barheight = 20)) +
#   scale_x_discrete(position = "top") +
#   theme_minimal() +
#   theme(text = element_text(size = 24),
#         axis.text.x = element_text(size = 28),
#         axis.title = element_blank(),
#         panel.grid = element_blank()) # 1000 by 1000
# 
# ggplot(factors_s2_blank2, aes(x = factor,
#                               y = reorder(capacity, desc(order)), fill = loading)) +
#   geom_tile(color = "black") +
#   # geom_text(aes(label = format(round(loading, 2), nsmall = 2))) +
#   scale_fill_distiller(palette = "RdYlBu", limits = c(-1, 1), breaks = c(-1, 0, 1),
#                        guide = guide_colorbar(title = element_blank(),
#                                               barheight = 20)) +
#   scale_x_discrete(position = "top") +
#   theme_minimal() +
#   theme(text = element_text(size = 24),
#         axis.text.x = element_text(size = 28),
#         axis.title = element_blank(),
#         panel.grid = element_blank()) # 1000 by 1000
# 
# ggplot(factors_s2_blank3, aes(x = factor,
#                               y = reorder(capacity, desc(order)), fill = loading)) +
#   geom_tile(color = "black") +
#   # geom_text(aes(label = format(round(loading, 2), nsmall = 2))) +
#   scale_fill_distiller(palette = "RdYlBu", limits = c(-1, 1), breaks = c(-1, 0, 1),
#                        guide = guide_colorbar(title = element_blank(),
#                                               barheight = 20)) +
#   scale_x_discrete(position = "top") +
#   theme_minimal() +
#   theme(text = element_text(size = 24),
#         axis.text.x = element_text(size = 28),
#         axis.title = element_blank(),
#         panel.grid = element_blank()) # 1000 by 1000
ggplot(factors_s2_long, aes(x = factor,
                            y = reorder(capacity, desc(order)), fill = loading)) +
  geom_tile(color = "black") +
  geom_text(aes(label = format(round(loading, 2), nsmall = 2)), size = 6) +
  scale_fill_distiller(palette = "RdYlBu", limits = c(-1, 1), breaks = c(-1, 0, 1),
                       guide = guide_colorbar(title = element_blank(),
                                              barheight = 20)) +
  scale_x_discrete(position = "top") +
  # geom_rect(aes(xmin = 0.51, xmax = 1.49, ymin = 14.55, ymax = 20.45),
  #           alpha = 0, color = "black", size = .5) +
  # geom_rect(aes(xmin = 1.51, xmax = 2.49, ymin = 6.55, ymax = 14.45),
  #           alpha = 0, color = "black", size = .5) +
  # geom_rect(aes(xmin = 2.51, xmax = 3.49, ymin = 0.55, ymax = 6.45),
  #           alpha = 0, color = "black", size = .5) +
  # theme_bw() +
  theme_minimal() +
  theme(text = element_text(size = 24),
        axis.text.x = element_text(size = 28),
        axis.title = element_blank(),
        panel.grid = element_blank()) # 1000 by 1000

Mean ratings

Studies 1-2

Mean ratings of 40 mental capacities for the 2 entities included in Studies 1-2. Participants responded on a 3-point scale (0 = “no”, 0.5 = “kinda”, 1 = “yes”). Error bars are nonparametric bootstrapped 95% confidence intervals. Mental capacities are grouped according to their dominant factor loading in Study 1 (adults).

# make dataframe
s12_plotting <- char_plotting %>%
  filter(study %in% c("study 1", "study 2")) %>%
  distinct()
# plot! (ordered by study 1 factor loadings)
s12 <- ggplot(s12_plotting,
               aes(y = Mean, x = reorder(wording, desc(s1_order)),
                   colour = factor(s1_color), shape = study)) +
  geom_point(stat = "identity", position = position_dodge(width = 0.6), size = 2) +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.4,
                position = position_dodge(width = 0.6)) +
  facet_wrap(~ character) +
  theme_bw() +
  scale_y_continuous(name = "\nMean rating",
                     limits = c(0, 1),
                     breaks = c(0, 0.5, 1),
                     labels = c("0\n(no)", "0.5\n(kinda)", "1\n(yes)")) +
  scale_shape_discrete(name = "Study:",
                       labels = c("Study 1: adults", "Study 2: 7-9y")) +
  # scale_colour_brewer(name = "Factor:",
  #                     type = "qual", palette = 6,
  #                     guide = FALSE) +
  scale_colour_manual(name = "Factor:",
                      values = c("#E41A1C", "#4DAF4A", "#377EB8"),
                      labels = c("BODY", "MIND", "HEART")) +
  coord_flip() +
  theme(text = element_text(size = 9),
      axis.title.y = element_blank(),
      axis.text.y = element_text(face = "italic",
                                 colour = palette_s1),
      panel.grid.minor = element_blank(),
      legend.position = "right")
s12

Mean factor scores

Studies 1-2

r1b <- brm(score ~ character * factor * age_group + (1 | subid) , tempC,
            family = "gaussian", cores = n_cores)

Compiling the C++ model
Start sampling

starting worker pid=1176 on localhost:11418 at 13:16:46.683
starting worker pid=1184 on localhost:11418 at 13:16:46.951
starting worker pid=1192 on localhost:11418 at 13:16:47.198
starting worker pid=1200 on localhost:11418 at 13:16:47.439

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).

Gradient evaluation took 0.012079 seconds
1000 transitions using 10 leapfrog steps per transition would take 120.79 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).

Gradient evaluation took 0.015809 seconds
1000 transitions using 10 leapfrog steps per transition would take 158.09 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).

Gradient evaluation took 0.022298 seconds
1000 transitions using 10 leapfrog steps per transition would take 222.98 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).

Gradient evaluation took 0.025095 seconds
1000 transitions using 10 leapfrog steps per transition would take 250.95 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)

# r1c <- brm(score ~ character * factor * age_group + (factor | subid) , tempC,
#             family = "gaussian", cores = n_cores)
# summary(r1c)

# plot
ggplot(scores_s12_plotting %>%
         ungroup() %>%
         mutate(factor = factor(factor,
                                # labels = c("Social-emotional",
                                #            "Bodily",
                                #            "Perceptual-cognitive")),
                                labels = c("Heart",
                                           "Body",
                                           "Mind")),
                age_group = factor(age_group,
                                  # levels = c("adults", "children_79"),
                                  # labels = c("adults", "children"))),
                                  levels = c("children_79", "adults"),
                                  labels = c("children", "adults"))),
       aes(x = age_group, y = Mean, color = character, shape = character)) +
  facet_wrap("factor", ncol = 3) +
  theme_bw() +
  theme(text = element_text(size = 28),
        legend.position = "bottom") +
  geom_point(size = 5, position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = Lower, ymax = Upper),
                width = 0, position = position_dodge(width = 0.4)) +
  scale_shape_manual(values = c(19, 15)) +
  labs(title = "Factor scores by age group",
       # subtitle = "Adults (Study 1) vs. children (Study 2)\n",
       x = "Age group",
       y = "Mean factor score") # 1000 by 500

Factor scores by age

Study 1

scores_s1_plotting <- d1 %>% 
  select(subid, age, character) %>% 
  distinct() %>% 
  mutate(subid = paste(character, subid, sep = "_")) %>%
  full_join(efa_d1_rotatedN$scores %>% 
              data.frame() %>% 
              rownames_to_column("subid")) %>%
  mutate(character = factor(character)) %>%
  rename(score_F1 = MR1, score_F2 = MR2, score_F3 = MR3) %>%
  filter(!is.na(score_F1), !is.na(score_F2), !is.na(score_F3), !is.na(age)) %>%
  gather(factor, score, starts_with("score_")) %>%
  mutate(factor = factor(factor))

ggplot(scores_s1_plotting %>%
         ungroup() %>%
         mutate(factor = factor(factor,
                                levels = c("score_F1", "score_F2", "score_F3"),
                                # labels = c("Social-emotional",
                                #            "Bodily",
                                #            "Perceptual-cognitive"))),
                                labels = c("Heart",
                                           "Body",
                                           "Mind"))),
       aes(x = age, y = score, color = character, fill = character, shape = character)) +
  facet_wrap("factor", ncol = 3) +
  theme_bw() +
  theme(text = element_text(size = 28),
        legend.position = "bottom") +
  # geom_smooth(method = "loess", alpha = 0.4) +
  geom_smooth(method = "lm", alpha = 0.4) +
  geom_point(size = 2) +
  scale_shape_manual(values = c(19, 15)) +
  labs(title = "Factor scores by adults' age",
       # subtitle = "Adults (Study 1)\n",
       x = "Age (years)",
       y = "Factor score") # 1000 by 500

Study 2

scores_s2_plotting <- d2 %>% 
  select(subid, age, character) %>% 
  distinct() %>% 
  mutate(subid = paste(character, subid, sep = "_")) %>%
  # full_join(efa_d2_rotatedN$scores %>% 
  full_join(efa_d12_all_rotated$scores %>% 
              data.frame() %>% 
              rownames_to_column("subid")) %>%
  mutate(character = factor(character)) %>%
  rename(score_F1 = MR1, score_F2 = MR2, score_F3 = MR3) %>%
  filter(!is.na(score_F1), !is.na(score_F2), !is.na(score_F3), !is.na(age)) %>%
  gather(factor, score, starts_with("score_")) %>%
  mutate(factor = factor(factor))

contrasts(scores_s2_plotting$factor) = cbind(factor1 = c(1, -1, 0),
                                factor3 = c(0, -1, 1))

contrasts(scores_s2_plotting$character) = cbind(robot = c(-1, 1))

r2 <- lmer(score ~ character * factor * scale(age, scale = F) + (1 | subid) , scores_s2_plotting)
summary(r2)

# r2b <- brm(score ~ character * factor * age + (1 | subid) , scores_s2_plotting %>% mutate(age = scale(age, scale = F)), family = "gaussian", cores = n_cores)
# summary(r2b)

# look at heart only
# r2c_heart <- brm(score ~ character * age, 
#            data = scores_s2_plotting %>% 
#              filter(factor == "score_F1") %>% 
#              mutate(age = scale(age, scale = F)), 
#            family = "gaussian", cores = n_cores)
# summary(r2c_heart)

# look at body only
# r2c_body <- brm(score ~ character * age, 
#            data = scores_s2_plotting %>% 
#              filter(factor == "score_F2") %>% 
#              mutate(age = scale(age, scale = F)), 
#            family = "gaussian", cores = n_cores)
# summary(r2c_body)

# look at mind only
# r2c_mind <- brm(score ~ character * age, 
#            data = scores_s2_plotting %>% 
#              filter(factor == "score_F3") %>% 
#              mutate(age = scale(age, scale = F)), 
#            family = "gaussian", cores = n_cores)
# summary(r2c_mind)

ggplot(scores_s2_plotting %>%
         ungroup() %>%
         mutate(factor = factor(factor,
                                labels = c("Heart",
                                           "Body",
                                           "Mind"))),
                         # labels = c("Social-emotional",
                         #            "Bodily",
                         #            "Perceptual-cognitive"))),
       aes(x = age, y = score, color = character, fill = character, shape = character)) +
  facet_wrap("factor", ncol = 3) +
  theme_bw() +
  theme(text = element_text(size = 28),
        legend.position = "bottom") +
  # geom_smooth(method = "loess", alpha = 0.4) +
  geom_smooth(method = "lm", alpha = 0.4) +
  geom_point(size = 2) +
  scale_shape_manual(values = c(19, 15)) +
  labs(title = "Factor scores by children's age",
       # subtitle = "Children (Study 2)\n",
       x = "Age (years)",
       y = "Factor score") # 1000 by 500

Studies 1-2

tempA <- d2 %>% 
  select(subid, age, character) %>% 
  distinct() %>% 
  mutate(subid = paste(character, subid, sep = "_")) %>%
  full_join(efa_d12_all_rotatedN$scores %>% 
              data.frame() %>% 
              rownames_to_column("subid")) %>%
  mutate(character = factor(character)) %>%
  rename(score_F1 = MR1, score_F2 = MR2, score_F3 = MR3) %>%
  filter(!is.na(score_F1), !is.na(score_F2), !is.na(score_F3), !is.na(age)) %>%
  gather(factor, score, starts_with("score_")) %>%
  mutate(factor = factor(factor,
                         labels = c("Heart",
                                    "Body",
                                    "Mind")))
                         # labels = c("Social-emotional",
                         #            "Bodily",
                         #            "Perceptual-cognitive")))

tempB <- scores_s12_plotting %>%
  filter(age_group == "adults") %>%
  ungroup() %>%
  mutate(factor = factor(factor,
                         labels = c("Heart",
                                    "Body",
                                    "Mind")),
                         # labels = c("Social-emotional",
                         #            "Bodily",
                         #            "Perceptual-cognitive")),
         age = 11) 

ggplot(tempA,
       aes(x = age, y = score, color = character, fill = character, shape = character)) +
  facet_wrap("factor", ncol = 3) +
  theme_bw() +
  theme(text = element_text(size = 28),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  geom_hline(data = tempB, aes(yintercept = Mean, color = character), lty = 2) +
  # geom_smooth(method = "loess", alpha = 0.4) +
  geom_smooth(method = "lm", alpha = 0.4) +
  geom_point(size = 2) +
  geom_point(data = tempB, aes(y = Mean),
             size = 4, position = position_dodge(width = 0.6)) +
  geom_errorbar(data = tempB, aes(ymin = Lower, ymax = Upper, y = Mean), width = 0,
                position = position_dodge(width = 0.6)) +
  scale_shape_manual(values = c(19, 15)) +
  scale_x_continuous(breaks = c(7:11), labels = c("7y", "8y", "9y", "10y", "adults")) +
  labs(title = "Factor scores by age",
       # subtitle = "Children (Study 2)\n",
       x = "Age",
       y = "Factor score") # 1000 by 500